---
title: "Mass valuation of agricultural land in Ukraine: regression results"
author: "Eduard Bukin ebukin@worldbank.org"
date: "November 20, 2024"
format: 
  docx:
    reference-doc: template.docx
editor: source
editor_options: 
  chunk_output_type: console
---

```{r echo=FALSE, warning=FALSE, message=FALSE}
library(here)
library(tidyverse)
library(sf)
library(arrow)
library(flextable)
library(glue)
# library(DescTools)
library(modelsummary)
library(scales)
library(ggExtra)
library(gtools)
library(broom)
# library(fixest)
library(rsample)
library(yardstick)
library(pins)

knitr::opts_chunk$set(
  warning = FALSE,
  error = FALSE, 
  message = FALSE, 
  echo = FALSE
  )

here("supporting-functions.R") |> source()

bd_paper <- board_folder(here("output", "paper-data"), versioned = TRUE)

FitFlextableToPage <- function(ft, pgwidth = 6.5, size = 9) {
  ft_out <-
    ft |> fontsize(size = size, part = "all") |> 
    font(
      fontname = "Times New Roman",
      part = "all") |> 
    autofit()
  
  ft_out <- width(ft_out,
                  width = dim(ft_out)$widths * pgwidth / (flextable_dim(ft_out)$widths))
  
  ft_out <- 
    ft_out |> 
    bold(part = "header") |> 
    valign(valign  = "center", part = "all") |> 
    align(j = -1, align = "center", part = "header") |> 
    align(j = -1, part = "body")
  
  return(ft_out)
}
```

\newpage

```{r table-1}
```

**Table 1. Descriptive statistics for transacted and non-transacted parcels**

```{r table-1-sum-stat}
ft_out <- 
  bd_paper |> 
  pin_read("dta-summary-stats") |> 
  mutate(across(where(is.character), ~ ifelse(is.na(.), " ", .))) |> 
  rename(
    `Transacted parcels` = reg_dta,
    `Not transacted parcels` = oos_dta,
    `Difference` = `Diff.`
  ) |> 
  flextable() |>
  add_footer_lines(
           value = as_paragraph(as_i("Note: "), "Statistics is reported as: Mean (SD), Column 3 reports the difference in means between parcels that were and were not transacted together with results of a two-sided Welch t-test assuming unequal variances for the significance of the difference. Significance levels are `*` < 0.05, '**' < 0.01, and '***' < 0.001")
  ) |> 
  FitFlextableToPage()

ft_out
```

\newpage

```{r table-2}
```

**Table 2. Regression results**

```{r table-2-reg-results}
reg_res_tbl <- bd_paper |> pin_read("tbl-all-mods")
reg_col_names <- reg_res_tbl$header$col_keys
# reg_col_names <-
#   ifelse(str_detect(reg_col_names, "4|5"),
#          str_c(reg_col_names, " FE"),
#          reg_col_names)
reg_res_tbl |>
  set_header_labels(values = reg_col_names) |> 
  FitFlextableToPage(size = 7) |> 
  width(j = 1, width = 1.5)
```

\newpage

```{r table-3}
```

**Table 3. Frequency of FE used for in-sample predictions**

```{r table-3-fe-in-sample-freq}
tbl3 <-
  pin_read(bd_paper, "tbl-fe-freq") |> 
  mutate(across(where(is.numeric), ~ ifelse(is.na(.), 0, .))) |> 
  mutate(`Total` = pick(`Cad. block`:Oblast) |> rowSums()) |> 
  filter(sample == "in") |> 
  select(-sample) |> 
  select(Region = region, Total, everything())
tbl3_country <- tbl3 |> filter(Region == "Country")

tbl3 |> filter(Region != "Country") |> 
  bind_rows(tbl3_country) |> 
  mutate(across(where(is.numeric), ~ number(., 1, big.mark = ","))) |> 
  bind_rows(
    tbl3_country |>
      mutate(
        Region = "Percentage",
        across(where(is.numeric), ~ . / tbl3_country$Total[[1]]),
        across(where(is.numeric), ~ percent(., .1)),
        Total = ""
      )
  ) |> 
  flextable() |> 
  FitFlextableToPage()
```


```{r gof-ensabmble}
reg_dta7 <- bd_paper |> pin_read("tbl-ensamble-r2") 

format_numbers <- function(dta) {
  dta |> 
    mutate(
      across(contains("r2"),  ~ number(., 0.001)),
      
      across(contains("nn"),  ~ number(., 1, big.mark = "")),
      across(matches("_p\\d{1,}"),  ~ number(., 0.01, big.mark = "")),
      across(where(is.numeric),  ~ number(., 0.1, big.mark = ""))
    )
}

tbl_actual_dta_sum <-
  reg_dta7 |> 
  select(id, contains("splits")) |>
  mutate(test_dta = map(splits, testing)) |> 
  bind_rows(
    tibble(id = "Full sample",
           test_dta = list(
             bind_rows(training(reg_dta7$splits[[1]]), 
                       testing(reg_dta7$splits[[1]]))
           ))
  ) |> 
  mutate(
    test_dta = map(
      test_dta,
      ~{
        .x |> 
          summarise(
            nn = n(),
            pr_p0 = exp(lprice_uah_2024_usd_ha_win) |> 
              stats::quantile(prob = 0, na.rm = TRUE),
            pr_p1 = exp(lprice_uah_2024_usd_ha_win) |> 
              stats::quantile(prob = 0.01, na.rm = TRUE),
            pr_p5 = exp(lprice_uah_2024_usd_ha_win) |> 
              stats::quantile(prob = 0.05, na.rm = TRUE),
            pr_p50 = exp(lprice_uah_2024_usd_ha_win) |> 
              stats::quantile(prob = 0.5, na.rm = TRUE),
            pr_p99 = exp(lprice_uah_2024_usd_ha_win) |> 
              stats::quantile(prob = 0.99, na.rm = TRUE),
            pr_p100 = exp(lprice_uah_2024_usd_ha_win) |> 
              stats::quantile(prob = 1, na.rm = TRUE),
            
            ac_p0 = exp(lprice_uah_2024_usd_ha) |> 
              stats::quantile(prob = 0, na.rm = TRUE),
            ac_p1 = exp(lprice_uah_2024_usd_ha) |> 
              stats::quantile(prob = 0.01, na.rm = TRUE),
            ac_p5 = exp(lprice_uah_2024_usd_ha) |> 
              stats::quantile(prob = 0.05, na.rm = TRUE),
            ac_p50 = exp(lprice_uah_2024_usd_ha) |> 
              stats::quantile(prob = 0.5, na.rm = TRUE),
            ac_p99 = exp(lprice_uah_2024_usd_ha) |> 
              stats::quantile(prob = 0.99, na.rm = TRUE),
            ac_p100 = exp(lprice_uah_2024_usd_ha) |> 
              stats::quantile(prob = 1, na.rm = TRUE),
            
            ar_p0 = area_ha |> stats::quantile(prob = 0, na.rm = TRUE),
            ar_p1 = area_ha |> stats::quantile(prob = 0.01, na.rm = TRUE),
            ar_p5 = area_ha |> stats::quantile(prob = 0.05, na.rm = TRUE),
            ar_p50 = area_ha |> stats::quantile(prob = 0.5, na.rm = TRUE),
            ar_p99 = area_ha |> stats::quantile(prob = 0.99, na.rm = TRUE),
            ar_p100 = area_ha |> stats::quantile(prob = 1, na.rm = TRUE)
          ) 
      }
    )
  ) |> 
  unnest(test_dta) |> select(-splits) |>  format_numbers() |> 
  pivot_longer(cols = c(everything(), -id)) |> 
  separate(name, into = c("name", "stats")) |> 
  left_join(tibble(name = c("nn", "pr", "ar", "ac"),
                   name_2 = c("Number of obs.", "Pred. price/ha", "Area, ha",
                              "Actual price/ha")),
            by = join_by(name)) |> 
  select(id, var = name_2, stats, value)


tbl_singe_mod_sum <-
  reg_dta7 |> 
  select(id, contains("predict_summ")) |> 
  unnest(predict_summ) |> 
  select(-contains("n_"), -nn) |> 
  format_numbers() |> 
  pivot_longer(c(everything(), -id, -model)) |> 
  separate(name, into = c("name", "stats")) |> 
  left_join(tibble(name = c("nn", "pr", "r2"),
                   name_2 = c("N obs.", "Pred. price/ha", "R2")),
            by = join_by(name))|> 
  select(id, model, var = name_2, stats, value)

tbl_averages_sum <-
  reg_dta7 |> 
  select(id, contains("predict_test_means")) |> 
  unnest(predict_test_means) |> 
  group_by(id) |> 
  summarise(
    r2 =  rsq_vec(lprice_uah_2024_usd_ha_win , p_price),
    pr_p0 = exp(p_price) |> stats::quantile(prob = 0, na.rm = TRUE),
    pr_p1 = exp(p_price) |> stats::quantile(prob = 0.01, na.rm = TRUE),
    pr_p5 = exp(p_price) |> stats::quantile(prob = 0.05, na.rm = TRUE),
    pr_p50 = exp(p_price) |> stats::quantile(prob = 0.5, na.rm = TRUE),
    pr_p99 = exp(p_price) |> stats::quantile(prob = 0.99, na.rm = TRUE),
    pr_p100 = exp(p_price) |> stats::quantile(prob = 1, na.rm = TRUE)
    ) |> 
  format_numbers() |> 
  mutate(model = "averages") |> 
  pivot_longer(c(everything(), -id, -model)) |> 
  separate(name, into = c("name", "stats")) |> 
  left_join(tibble(name = c("nn", "pr", "r2"),
                   name_2 = c("N obs.", "Pred. price/ha", "R2")),
            by = join_by(name)) |> 
  select(id, model, var = name_2, stats, value)

ensamble_sum <-
  reg_dta7 |> 
  select(id, contains("ensamble")) |>
  select(id, contains("_summ")) |>
  pivot_longer(contains("_summ"), names_to = "model") |> 
  unnest(value) |> 
  mutate(across(contains("n_adm"), ~ percent(. / nn, 0.1))) |> 
  format_numbers() |> 
  select(-nn) |> 
  pivot_longer(c(everything(), -id, -model))  |> 
  separate(name, into = c("name", "stats")) |> 
  left_join(tibble(name = c("pr", "r2", "N"),
                   var = c("Pred. price/ha", "R2", "% of FE from")),
            by = join_by(name)) |> 
  select(-name)
  
tbl_sum_all <- 
  tbl_actual_dta_sum |> mutate(type = "Testing data") |> 
  bind_rows(
    tbl_singe_mod_sum |> 
      left_join(
        tibble(model = c("adm1_id", "adm2_id", "adm2old_id", "adm3_id", "adm6_id", "adm7_id", "adm8_id"),
               type = c("Oblast", "Rayon", "Old rayon", "Community",
                        "Village", "Cad. zone", "Cad. block"))
      ) |> 
      select(-model)
  ) |> 
  bind_rows(
    ensamble_sum |> 
      left_join(
        tibble(model = c("ensamble_5_summ", "ensamble_10_summ", "ensamble_40_summ"),
               type = c("Ensamble 5", "Ensamble 10", "Ensamble 40"))
      ) |> 
      filter(!is.na(type)) |> 
      select(-model)
  )  |> 
  bind_rows(
    tbl_averages_sum |> 
      left_join(tibble(model = c("averages"), type = c("Reg. averages"))) |> 
      select(-model)
  )  |> 
  pivot_wider(names_from = id, values_from = value)


tbl_sum_all <- 
  tbl_sum_all |> 
  left_join(
    tibble(stats = c(NA, "p0", "p1", "p5", "p50", "p99", "p100"),
           stats2 = c(NA, "Min", "Q1", "Q5", "Median", "Q99", "Max")),
    by = join_by(stats)
    ) |> 
  mutate(stats = ifelse(is.na(stats2), stats, stats2)) |> 
  select(-stats2)
```


\


```{r table-4}
```

**Table 4. Predictive power of different models by cross-validation split**

```{r table-4-cross-valid-predictive-power}
tbl_sum_all |> filter(is.na(stats)) |> filter(var == "R2") |> 
  select(Model = type, var, everything()) |> select(-stats, -var) |> 
  select(-"Full sample") |> 
  flextable() |> FitFlextableToPage()
```

\

```{r table-5}
```

**Table 5: Prediction table by region**

```{r table-5-outsample-country}
tbl5 <- 
  bd_paper |> pin_read("tbl-region-stat") |> 
  rename(name2 = name) |> 
  pivot_longer(c(everything(), -name2)) |> 
  pivot_wider(names_from = name2, values_from = value) |> 
  mutate(
    actual_nmv = p_actual / p_nmv - 1,
    reg_nmv = p_reg  / p_nmv - 1,
    aver_nmv = p_kvaverage  / p_nmv - 1
  ) |> 
  mutate(
    across(c(n , n_nmv, contains("p_")), ~number(., 1, big.mark = ",")),
    across(c(area, size), ~number(., .01, big.mark = ",")), 
    across(where(is.numeric), ~number(. * 100, 0.1, big.mark = ",")),
    across(where(is.character), ~ifelse(is.na(.), "", .))
  ) |> 
  rename(name2 = name) |> 
  pivot_longer(c(everything(), -name2)) |> 
  pivot_wider(names_from = name2, values_from = value) |> 
  select(name, Total, Sample = `In-sample`, Prediction = `Out-of-sample`,
         everything())

tribble(
  ~name, ~` `,
  "n", "No. of parcels ",
  "n_nmv", "  w. NMV",
  "area",  "Area, 1000 ha",
  "size", "Plot size, ha",
  "p_actual",  "Actual price",
  "p_reg", "Pred. price",
  "p_kvaverage", "Region avg",
  "p_nmv", "NMV",
  "aver_nmv", "Reg. avg/NMV",
  "reg_nmv", "Pred/NMV",
  "actual_nmv",  "Actual/NMV"
) |> 
  
  left_join(tbl5, by = join_by(name)) |> 
  select(-name) |> 
  flextable() |>
  FitFlextableToPage() |>
  width(j = 1, width = 1.25)
```

\newpage

```{r figure-1}
```


![](output/paper-img/gof-model-restricted.png){}

**Figure 1. Actual vs. predicted prices by model**

\

```{r figure-2}
```

\newpage

**Appendix figures and tables**

\

```{r figure-A1}
```

![](output/paper-img/plots-size-price.png)

\

**Figure A1. Parcel price per ha versus size relationship**

\

```{r figure-A2}
```

![](output/paper-img/distribution-area-nmv-soil.png){}

\

**Figure A2. Plot area, NMV and soil characteristics distribution in regression and out-of-sample data**

\

```{r figure-A3}
```

![](output/paper-img/distribution-distances.png){}

\

**Figure A3. Distances to infrastructure distribution in regression and out-of-sample data**

\


```{r figure-A4-kvartal}
```

![](output/paper-img/trans-freq-adm8_id.png){}

\

**Figure A4. Transactions frequency by cadastral kvartal**

\



```{r figure-A5-zone}
```

![](output/paper-img/trans-freq-adm7_id.png){}

\

**Figure A5. Transactions frequency by cadastral zone**

\



```{r figure-A6-village}
```

![](output/paper-img/trans-freq-adm6_id.png){}

\

**Figure A6. Transactions frequency by village**

\



```{r figure-A7-hromada}
```

![](output/paper-img/trans-freq-adm3_id.png){}

\

**Figure A7. Transactions frequency by hromada**

\



```{r figure-A8-old-rayon}
```

![](output/paper-img/trans-freq-adm2old_id.png){}

\

**Figure A8. Transactions frequency by old rayon**

\



```{r figure-A9-rayon}
```

![](output/paper-img/trans-freq-adm2_id.png){}

\

**Figure A9. Transactions frequency by rayon**

\



```{r figure-A10-nmv-freq}
```

![](output/paper-img/nmv-freq-hromada.png){}

\


**Figure A10. Frequency of NMV observations by hromada**

\


\newpage

```{r table-A1}
```

**Appendix table A1: Sizes and prices outliers in regression data and by cross validation split**

```{r table-A1-outliers}
tbl_sum_all |> 
  filter(!is.na(stats), !stats %in% c("adm1")) |> 
  filter(str_detect(var, "Area|Actu|Pred"), type == "Testing data") |>
  mutate(
    var = ifelse(var == "Pred. price/ha", "Actual price (winsorized)/ha", var) |> 
      str_replace("/ha", " 2024 USD/ha")) |>
  select(-type) |> rename(`  ` = stats) |>  
  as_grouped_data(groups = "var") |> 
  as_flextable(hide_grouplabel = TRUE) |>
  bold(j = 1, i = ~ !is.na(var), bold = TRUE, part = "body") |> 
  FitFlextableToPage()
```

\newpage

```{r table-A2}
```

**Appendix table A2: Covariates contribution to variance explained in percent**

```{r table-3-r2-contrib}
pin_read(bd_paper, "tbl-r2-contrib")$full_model |> 
  rename(`  ` = var, `N variables` = n_var) |> 
  rename_with(~str_remove_all(., "\\(\\d\\) ")) |> 
  rename_with(~str_remove_all(., " FE")) |> 
  mutate(across(where(is.numeric), ~percent(., 0.1))) |> 
  flextable() |> FitFlextableToPage() |>
  width(j = 1, width = 1.8)
```

\newpage

```{r table-A3}
```

**Appendix table A3: Frequencies of the FE used by region in out-of-sample prediction**

```{r table-3-fe-oos-freq}
tbl3a <-
  pin_read(bd_paper, "tbl-fe-freq") |> 
  mutate(across(where(is.numeric), ~ ifelse(is.na(.), 0, .))) |> 
  mutate(`Total` = pick(`Cad. block`:Oblast) |> rowSums()) |> 
  filter(sample == "out") |> 
  select(-sample) |> 
  select(Region = region, Total, everything())
tbl3a_country <- tbl3a |> filter(Region == "Country")

tbl3a |> filter(Region != "Country") |> 
  bind_rows(tbl3a_country) |> 
  mutate(across(where(is.numeric), ~ number(., 1, big.mark = ","))) |> 
  bind_rows(
    tbl3a_country |>
      mutate(
        Region = "Percentage",
        across(where(is.numeric), ~ . / tbl3a_country$Total[[1]]),
        across(where(is.numeric), ~ percent(., .1)),
        Total = ""
      )
  ) |> 
  flextable() |> 
  FitFlextableToPage()
```
